
rm(list=ls())
library(tidyverse)
library(scales)
library(broom)
#dt<-import data12.csv as dt


# Panel A ----
fit_c0 <- lm( # A.
  c_eval ~ c_exp + p_exp +
    female + educ + pidc + pid_dk + age + black + hispanic + asian,
  weights = weight, data = dt
)

panel_a <- tidy(fit_c0) |>
  filter(term %in% c("c_exp","p_exp")) |>
  dplyr::select(term, estimate, std.error) |>
  mutate(
    ci_low = estimate - 1.96*std.error,
    ci_up  = estimate + 1.96*std.error,
    Evaluations = "Court Evaluations",
    terms = ifelse(term == "p_exp", "Police", "Courts")
  )


# panel B ----

fit_c1int <- lm( # B.
  p_eval ~ c_exp + p_exp  +cnp*c_exp+
    female + educ + pidc + pid_dk + age + black + hispanic + asian,
  weights = weight, data = dt
)

panel_b <- tidy(fit_c1int) |>
  filter(term %in% c("p_exp","c_exp")) |>
  dplyr::select(term, estimate, std.error) |>
  mutate(
    ci_low = estimate - 1.96*std.error,
    ci_up  = estimate + 1.96*std.error,
    Evaluations = "Police Evaluations (conditional court effect)",
    terms = ifelse(term == "p_exp", "Police", "Courts"))


b  <- coef(fit_c1int)
V  <- vcov(fit_c1int)

# names: "p_exp", "pnc", "p_exp:pnc"
Delta_hat <- b["c_exp"] + b["cnp"] + b["c_exp:cnp"]

Var_Delta <- V["c_exp","c_exp"] +
  V["cnp","cnp"] +
  V["c_exp:cnp","c_exp:cnp"] +
  2*( V["c_exp","cnp"] + V["c_exp","c_exp:cnp"] + V["cnp","c_exp:cnp"] )

SE_Delta  <- sqrt(Var_Delta)
CI        <- c(Delta_hat - 1.96*SE_Delta, Delta_hat + 1.96*SE_Delta)

md<-data.frame(term="cnc",estimate = Delta_hat, se = SE_Delta, ci_low = CI[1], ci_high = CI[2])
colnames(md)<-c("term","estimate","std.error","ci_low","ci_up")
md$Evaluations<-c("Police Evaluations (conditional court effect)");md$terms<-c("Courts")

panel_b<-rbind(panel_b,md)




# panel C ----
fit_p0 <- lm( # C.
  p_eval ~ c_exp + p_exp +
    female + educ + pidc + pid_dk + age + black + hispanic + asian,
  weights = weight, data = dt
)

panel_c <- tidy(fit_p0) |>
  filter(term %in% c("c_exp","p_exp")) |>
  dplyr::select(term, estimate, std.error) |>
  mutate(
    ci_low = estimate - 1.96*std.error,
    ci_up  = estimate + 1.96*std.error,
    Evaluations = "Police Evaluations",
    terms = ifelse(term == "p_exp", "Police", "Courts")
  )

# panel D ----
fit_p1int <- lm( # D.
  c_eval ~ c_exp + p_exp  +pnc*p_exp+
    female + educ + pidc + pid_dk + age + black + hispanic + asian,
  weights = weight, data = dt
)

panel_d <- tidy(fit_p1int) |>
  filter(term %in% c("c_exp","p_exp")) |>
  dplyr::select(term, estimate, std.error) |>
  mutate(
    ci_low = estimate - 1.96*std.error,
    ci_up  = estimate + 1.96*std.error,
    Evaluations = "Court Evaluations (conditional police effect)",
    terms = ifelse(term == "p_exp", "Police", "Courts")
  )

b  <- coef(fit_p1int)
V  <- vcov(fit_p1int)

Delta_hat <- b["p_exp"] + b["pnc"] + b["p_exp:pnc"]

Var_Delta <- V["p_exp","p_exp"] +
  V["pnc","pnc"] +
  V["p_exp:pnc","p_exp:pnc"] +
  2*( V["p_exp","pnc"] + V["p_exp","p_exp:pnc"] + V["pnc","p_exp:pnc"] )

SE_Delta  <- sqrt(Var_Delta)
CI        <- c(Delta_hat - 1.96*SE_Delta, Delta_hat + 1.96*SE_Delta)

md<-data.frame(term="pnc",estimate = Delta_hat, se = SE_Delta, ci_low = CI[1], ci_high = CI[2])
colnames(md)<-c("term","estimate","std.error","ci_low","ci_up")
md$Evaluations<-c("Court Evaluations (conditional police effect)");md$terms<-c("Police")

panel_d<-rbind(panel_d,md)


# Plot figure 1 ----

plotdf <- dplyr::bind_rows(panel_a, panel_b,panel_c, panel_d) %>%
  dplyr::mutate(
    Evaluations = factor(
      Evaluations,
      levels = c(
        "Court Evaluations",
        "Court Evaluations (conditional police effect)",
        "Police Evaluations (conditional court effect)",
        "Police Evaluations"
      )
    ),
    terms = factor(terms, levels = c("Courts","Police"))
  )

xlim_pad <- range(plotdf$ci_low, plotdf$ci_up, na.rm = TRUE)
xlim_pad <- xlim_pad + c(-0.02, 0.02)

mid_lab <- "Court Evaluations (conditional police effect)"
mid_lab1 <- "Police Evaluations (conditional court effect)"
plotdf <- plotdf %>%
  group_by(Evaluations, terms) %>%
  mutate(
    terms_new = case_when(
      Evaluations == mid_lab & terms == "Police" & row_number() == 1 ~ "Police",
      Evaluations == mid_lab & terms == "Police" & row_number() == 2 ~ "Police only",
      Evaluations == mid_lab1 & terms == "Courts" & row_number() == 1 ~ "Courts",
      Evaluations == mid_lab1 & terms == "Courts" & row_number() == 2 ~ "Courts only",
      TRUE ~ terms
    ),
    y_off = case_when(
      Evaluations == mid_lab & terms == "Police" & row_number() == 1 ~ -0.10,
      Evaluations == mid_lab & terms == "Police" & row_number() == 2 ~  0.10,
      Evaluations == mid_lab1 & terms == "Courts" & row_number() == 1 ~ -0.10,
      Evaluations == mid_lab1 & terms == "Courts" & row_number() == 2 ~  0.10,
      TRUE ~ 0
    )
  ) %>% ungroup()

plotdf$Evaluations<-ifelse(plotdf$Evaluations=="Court Evaluations","A. Court Evaluations",
                    ifelse(plotdf$Evaluations=="Court Evaluations (conditional police effect)","B. Court Evaluations (conditional police effect)", 
                    ifelse(plotdf$Evaluations=="Police Evaluations","C. Police Evaluations","D. Police Evaluations (conditional court effect)"       
                           )))

plotdf$Evaluations<-factor(plotdf$Evaluations,levels=c("A. Court Evaluations","B. Court Evaluations (conditional police effect)","C. Police Evaluations","D. Police Evaluations (conditional court effect)"))

ggplot(
  plotdf,
  aes(y = terms, x = estimate, xmin = ci_low, xmax = ci_up,
      shape = terms_new, linetype = terms_new, color = terms_new)
) +
  geom_point(size = 3, position = position_nudge(y = plotdf$y_off)) +
  geom_errorbarh(height = .14, position = position_nudge(y = plotdf$y_off)) +
  geom_vline(xintercept = 0, linetype = 5, col = "grey20") +
  xlab("") + ylab("") +
  scale_x_continuous(limits = xlim_pad, n.breaks = 8, minor_breaks = 0) +
  # --- define values for ALL used levels, including "Police"
  scale_color_manual(
    values = c("Courts"="black","Police"="black","Police only"="black",
               "Police"="black","Police only"="black","Courts only"="black"),
    breaks = c("Courts","Police","Police only","Courts only"),   # single Police entry in legend
    labels = c("Courts","Police","Police only","Courts only")
  ) +
  scale_shape_manual(
    values = c("Courts"=0,"Police"=1,"Police only"=17,"Courts only"=18),
    breaks = c("Courts","Police","Police only","Courts only"),
    labels = c("Courts","Police","Police only","Courts only")
  ) +
  scale_linetype_manual(
    values = c("Courts"="longdash","Police"="solid"
               ,"Police only"="dotted","Courts only"="dotted"),
    breaks = c("Courts","Police","Police only","Courts only"),
    labels = c("Courts","Police","Police only","Courts only")
  ) +
  guides(color = guide_legend(title = "Experiences"),
         shape = guide_legend(title = "Experiences"),
         linetype = guide_legend(title = "Experiences")) +
  facet_wrap(vars(Evaluations), nrow = 2, drop = FALSE) +
  theme_bw() +
  theme(legend.position = "bottom")









